﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using SpiceNet.Circuits;
using SpiceNet.Diagnostics;
using System.Numerics;
using SpiceNet.Simulations;
using SpiceNet.Integration;
using SpiceNet.Parameters;

namespace SpiceNet.Models
{
    /// <summary>
    /// A class for a Jfet instance
    /// </summary>
    public class Jfet : Instance
    {
        /// <summary>
        /// Get parameters
        /// </summary>
        public override Dictionary<string, Parameterized.IP> ParamTable
        {
            get { return _pTable; }
        }

        #region Parameters
        private enum _c
        {
            /* device parameters */
            JFET_AREA = 1,
            JFET_IC_VDS = 2,
            JFET_IC_VGS = 3,
            JFET_IC = 4,
            JFET_OFF = 5,
            JFET_TEMP = 6,
            /* device questions */
            JFET_DRAINNODE = 301,
            JFET_GATENODE = 302,
            JFET_SOURCENODE = 303,
            JFET_DRAINPRIMENODE = 304,
            JFET_SOURCEPRIMENODE = 305,
            JFET_VGS = 306,
            JFET_VGD = 307,
            JFET_CG = 308,
            JFET_CD = 309,
            JFET_CGD = 310,
            JFET_GM = 311,
            JFET_GDS = 312,
            JFET_GGS = 313,
            JFET_GGD = 314,
            JFET_QGS = 315,
            JFET_CQGS = 316,
            JFET_QGD = 317,
            JFET_CQGD = 318,
            JFET_CS = 319,
            JFET_POWER = 320
        }
        private static Dictionary<string, IP> _pTable = new Dictionary<string, IP>
        {
            { "off", new IP(IF.IOPU, (int)_c.JFET_OFF, T.FLAG, "Device initially off") },
            { "ic", new IP(IF.IOPAU, (int)_c.JFET_IC, T.REALVEC, "Initial VDS,VGS vector") },
            { "area", new IP(IF.IOPU, (int)_c.JFET_AREA, T.REAL, "Area factor") },
            { "ic-vds", new IP(IF.IOPAU, (int)_c.JFET_IC_VDS, T.REAL, "Initial D-S voltage") },
            { "ic-vgs", new IP(IF.IOPAU, (int)_c.JFET_IC_VGS, T.REAL, "Initial G-S volrage") },
            { "temp", new IP(IF.IOPU, (int)_c.JFET_TEMP, T.REAL, "Instance temperature") },
            { "drain-node", new IP(IF.OPU, (int)_c.JFET_DRAINNODE, T.INTEGER, "Number of drain node") },
            { "gate-node", new IP(IF.OPU, (int)_c.JFET_GATENODE, T.INTEGER, "Number of gate node") },
            { "source-node", new IP(IF.OPU, (int)_c.JFET_SOURCENODE, T.INTEGER, "Number of source node") },
            { "drain-prime-node", new IP(IF.OPU, (int)_c.JFET_DRAINPRIMENODE, T.INTEGER, "Internal drain node") },
            { "source-prime-node", new IP(IF.OPU, (int)_c.JFET_SOURCEPRIMENODE, T.INTEGER, "Internal source node") },
            { "vgs", new IP(IF.OP, (int)_c.JFET_VGS, T.REAL, "Voltage G-S") },
            { "vgd", new IP(IF.OP, (int)_c.JFET_VGD, T.REAL, "Voltage G-D") },
            { "ig", new IP(IF.OP, (int)_c.JFET_CG, T.REAL, "Current at gate node") },
            { "id", new IP(IF.OP, (int)_c.JFET_CD, T.REAL, "Current at drain node") },
            { "is", new IP(IF.OP, (int)_c.JFET_CS, T.REAL, "Source current") },
            { "igd", new IP(IF.OP, (int)_c.JFET_CGD, T.REAL, "Current G-D") },
            { "gm", new IP(IF.OP, (int)_c.JFET_GM, T.REAL, "Transconductance") },
            { "gds", new IP(IF.OP, (int)_c.JFET_GDS, T.REAL, "Conductance D-S") },
            { "ggs", new IP(IF.OP, (int)_c.JFET_GGS, T.REAL, "Conductance G-S") },
            { "ggd", new IP(IF.OP, (int)_c.JFET_GGD, T.REAL, "Conductance G-D") },
            { "qgs", new IP(IF.OPU, (int)_c.JFET_QGS, T.REAL, "Charge storage G-S junction") },
            { "qgd", new IP(IF.OPU, (int)_c.JFET_QGD, T.REAL, "Charge storage G-D junction") },
            { "cqgs", new IP(IF.OPU, (int)_c.JFET_CQGS, T.REAL, "Capacitance due to charge storage G-S junction") },
            { "cqgd", new IP(IF.OPU, (int)_c.JFET_CQGD, T.REAL, "Capacitance due to charge storage G-D junction") },
            { "p", new IP(IF.OPU, (int)_c.JFET_POWER, T.REAL, "Power dissipated by the JFET") }
        };

        /// <summary>
        /// Instance optional parameters
        /// </summary>
        public bool JFEToff { get; set; }
        public Parameter<double> JFETarea { get; private set; }
        public Parameter<double> JFETicVDS { get; private set; }
        public Parameter<double> JFETicVGS { get; private set; }
        public Parameter<double> JFETtemp { get; private set; }
        #endregion

        #region Constants
        private const int JFETvgs = 0;
        private const int JFETvgd = 1;
        private const int JFETcg = 2;
        private const int JFETcd = 3;
        private const int JFETcgd = 4;
        private const int JFETgm = 5;
        private const int JFETgds = 6;
        private const int JFETggs = 7;
        private const int JFETggd = 8;
        private const int JFETqgs = 9;
        private const int JFETcqgs = 10;
        private const int JFETqgd = 11;
        private const int JFETcqgd = 12;

        private const int JFETRDNOIZ = 0;
        private const int JFETRSNOIZ = 1;
        private const int JFETIDNOIZ = 2;
        private const int JFETFLNOIZ = 3;
        private const int JFETTOTNOIZ = 4;
        #endregion

        #region Calculated constants
        public double JFETtSatCur { get; private set; } /* temperature adjusted saturation current */
        public double JFETtGatePot { get; private set; }    /* temperature adjusted gate potential */
        public double JFETtCGS { get; private set; }    /* temperature corrected G-S capacitance */
        public double JFETtCGD { get; private set; }    /* temperature corrected G-D capacitance */
        public double JFETcorDepCap { get; private set; }   /* joining point of the fwd bias dep. cap eq.s */
        public double JFETvcrit { get; private set; }   /* critical voltage for the instance */
        public double JFETf1 { get; private set; }      /* coefficient of capacitance polynomial exp */
        #endregion

        /// <summary>
        /// Private variables
        /// </summary>
        private int JFETstate; /* poprivate inter to start of state vector for jfet */
        private int JFETdrainNode;  /* number of drain node of jfet */
        private int JFETgateNode;   /* number of gate node of jfet */
        private int JFETsourceNode; /* number of source node of jfet */
        private int JFETdrainPrimeNode; /* number of private internal drain node of jfet */
        private int JFETsourcePrimeNode;    /* number of private internal source node of jfet */
        // private int JFETmode;
        private Noise.InstanceNoiseData[] JFETnVar = new Noise.InstanceNoiseData[5];

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        public Jfet(string name)
            : base(name, 3, typeof(JfetModel))
        {
            JFEToff = false;
            JFETarea = new Parameter<double>(1.0);
            JFETicVDS = new Parameter<double>();
            JFETicVGS = new Parameter<double>();
            JFETtemp = new Parameter<double>();
        }

        /// <summary>
        /// Set parameter
        /// </summary>
        protected override void Param(int id, object value, Circuit ckt = null)
        {
            double[] v = null;
            switch ((_c)id)
            {
                case _c.JFET_TEMP:
                    JFETtemp.Par((double)value + Circuit.CONSTCtoK);
                    break;
                case _c.JFET_AREA:
                    JFETarea.Par((double)value);
                    break;
                case _c.JFET_IC_VDS:
                    JFETicVDS.Par((double)value);
                    break;
                case _c.JFET_IC_VGS:
                    JFETicVGS.Par((double)value);
                    break;
                case _c.JFET_OFF:
                    JFEToff = (bool)value;
                    break;
                case _c.JFET_IC:
                    v = (double[])value;
                    switch (v.Length)
                    {
                        case 2:
                            JFETicVGS.Par(v[1]); goto case 1;
                        case 1:
                            JFETicVDS.Par(v[0]); break;
                        default:
                            throw new BadParameterException();
                    }
                    break;
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Ask a parameter
        /// </summary>
        /// <param name="id"></param>
        /// <param name="ckt"></param>
        /// <returns></returns>
        protected override object Ask(int id, Circuit ckt = null)
        {
            switch ((_c)id)
            {
                case _c.JFET_TEMP:
                    return JFETtemp - Circuit.CONSTCtoK;
                case _c.JFET_AREA:
                    return JFETarea.Value;
                case _c.JFET_IC_VDS:
                    return JFETicVDS.Value;
                case _c.JFET_IC_VGS:
                    return JFETicVGS.Value;
                case _c.JFET_OFF:
                    return JFEToff;
                case _c.JFET_DRAINNODE:
                    return JFETdrainNode;
                case _c.JFET_GATENODE:
                    return JFETgateNode;
                case _c.JFET_SOURCENODE:
                    return JFETsourceNode;
                case _c.JFET_DRAINPRIMENODE:
                    return JFETdrainPrimeNode;
                case _c.JFET_SOURCEPRIMENODE:
                    return JFETsourcePrimeNode;
                case _c.JFET_VGS:
                    return ckt.State.States[0][JFETstate + JFETvgs];
                case _c.JFET_VGD:
                    return ckt.State.States[0][JFETstate + JFETvgd];
                case _c.JFET_CG:
                    return ckt.State.States[0][JFETstate + JFETcg];
                case _c.JFET_CD:
                    return ckt.State.States[0][JFETstate + JFETcd];
                case _c.JFET_CGD:
                    return ckt.State.States[0][JFETstate + JFETcgd];
                case _c.JFET_GM:
                    return ckt.State.States[0][JFETstate + JFETgm];
                case _c.JFET_GDS:
                    return ckt.State.States[0][JFETstate + JFETgds];
                case _c.JFET_GGS:
                    return ckt.State.States[0][JFETstate + JFETggs];
                case _c.JFET_GGD:
                    return ckt.State.States[0][JFETstate + JFETggd];
                case _c.JFET_QGS:
                    return ckt.State.States[0][JFETstate + JFETqgs];
                case _c.JFET_CQGS:
                    return ckt.State.States[0][JFETstate + JFETcqgs];
                case _c.JFET_QGD:
                    return ckt.State.States[0][JFETstate + JFETqgd];
                case _c.JFET_CQGD:
                    return ckt.State.States[0][JFETstate + JFETcqgd];
                case _c.JFET_CS:
                    if (ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingAc))
                    {
                        throw new AcCurrentPowerException();
                    }
                    else
                    {
                        double value = -ckt.State.States[0][JFETstate + JFETcd];
                        value -= ckt.State.States[0][JFETstate + JFETcg];
                        return value;
                    }
                case _c.JFET_POWER:
                    if (ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingAc))
                    {
                        throw new AcCurrentPowerException();
                    }
                    else
                    {
                        double value = ckt.State.States[0][JFETstate + JFETcd] *
                                ckt.State.Solution[JFETdrainNode];
                        value += ckt.State.States[0][JFETstate + JFETcg] *
                                ckt.State.Solution[JFETgateNode];
                        value -= (ckt.State.States[0][JFETstate + JFETcd] +
                                ckt.State.States[0][JFETstate + JFETcg]) *
                        ckt.State.Solution[JFETsourceNode];
                        return value;
                    }
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Setup the Jfet
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void Setup(Circuit ckt, Model mod, object p = null)
        {
            JfetModel model = mod as JfetModel;
            List<NodeType> extra = new List<NodeType>();
            if (model.JFETsourceResist != 0.0)
                extra.Add(NodeType.Voltage);
            if (model.JFETdrainResist != 0.0)
                extra.Add(NodeType.Voltage);

            // Bind the nodes
            CircuitNode[] nodes = BindNodes(ckt, extra.ToArray());
            JFETdrainNode = nodes[0].Number;
            JFETgateNode = nodes[1].Number;
            JFETsourceNode = nodes[2].Number;
            int index = 3;
            JFETdrainPrimeNode = model.JFETdrainResist != 0.0 ? nodes[index++].Number : JFETdrainNode;
            JFETsourcePrimeNode = model.JFETsourceResist != 0.0 ? nodes[index++].Number : JFETsourceNode;

            // Create the states
            JFETstate = CreateState(ckt, 13);
        }

        /// <summary>
        /// Do temperature dependent calculations
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void Temperature(Circuit ckt, Model mod, object p = null)
        {
            JfetModel model = mod as JfetModel;
            double vt;
            double kt;
            double arg;
            double fact2;
            double egfet;
            double pbfact;
            double gmanew;
            double ratio1;
            double cjfact1;
            var tp = p as JfetModel.JfetTempParameters;

            if (!JFETtemp.Given)
                JFETtemp.Value = ckt.Config.Temperature;
            vt = JFETtemp * Circuit.CONSTKoverQ;
            fact2 = JFETtemp / Circuit.CONSTRefTemp;
            ratio1 = JFETtemp / model.JFETtnom - 1;
            JFETtSatCur = model.JFETgateSatCurrent * Math.Exp(ratio1 * 1.11 / vt);
            JFETtCGS = model.JFETcapGS * tp.cjfact;
            JFETtCGD = model.JFETcapGD * tp.cjfact;
            kt = Circuit.CONSTBoltz * JFETtemp;
            egfet = 1.16 - (7.02e-4 * JFETtemp * JFETtemp) /
                    (JFETtemp + 1108);
            arg = -egfet / (kt + kt) + 1.1150877 / (Circuit.CONSTBoltz * (Circuit.CONSTRefTemp + Circuit.CONSTRefTemp));
            pbfact = -2 * vt * (1.5 * Math.Log(fact2) + Circuit.CHARGE * arg);
            JFETtGatePot = fact2 * tp.pbo + pbfact;
            gmanew = (JFETtGatePot - tp.pbo) / tp.pbo;
            cjfact1 = 1 + .5 * (4e-4 * (JFETtemp - Circuit.CONSTRefTemp) - gmanew);
            JFETtCGS *= cjfact1;
            JFETtCGD *= cjfact1;

            JFETcorDepCap = model.JFETdepletionCapCoeff *
                    JFETtGatePot;
            JFETf1 = JFETtGatePot * (1 - Math.Exp((1 - .5) * tp.xfc)) / (1 - .5);
            JFETvcrit = vt * Math.Log(vt / (Circuit.CONSTroot2 * JFETtSatCur));
        }

        /// <summary>
        /// Load the JFet
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void Load(Circuit ckt, Model mod, object p = null)
        {
            JfetModel model = mod as JfetModel;

            double beta;
            double betap;
            double capgd;
            double capgs;
            double cd;
            double cdhat = 0.0;
            double cdrain;
            double cdreq;
            // double ceq;
            double ceqgd;
            double ceqgs;
            double cg;
            double cgd;
            double cghat = 0.0;
            double csat;
            double czgd;
            double czgdf2;
            double czgs;
            double czgsf2;
            double delvds;
            double delvgd;
            double delvgs;
            double evgd;
            double evgs;
            double fcpb2;
            double gdpr;
            double gds;
            // double geq;
            double ggd;
            double ggs;
            double gm;
            double gspr;
            double sarg;
            // double twob;
            double twop;
            double vds;
            double vgd;
            double vgdt;
            double vgs;
            double vgst;
            // double xfact;
            /* Modification for Sydney University JFET model */
            double vto; // phib
            double apart, cpart; // bpart
            double Bfac;
            // double vgd1, vgs1;
            /* end Sydney University mod. */
            bool icheck;
            bool ichk1;
            // int error;
            CircuitState state = ckt.State;

            /*
             *  dc model parameters 
             */
            beta = model.JFETbeta * JFETarea;
            gdpr = model.JFETdrainConduct * JFETarea;
            gspr = model.JFETsourceConduct * JFETarea;
            csat = JFETtSatCur * JFETarea;
            /*
             *    initialization
             */
            icheck = true;
            if ((ckt.Mode & Circuit.Modes.InitSmSig) != 0)
            {
                vgs = state.States[0][JFETstate + JFETvgs];
                vgd = state.States[0][JFETstate + JFETvgd];
            }
            else if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
            {
                vgs = state.States[1][JFETstate + JFETvgs];
                vgd = state.States[1][JFETstate + JFETvgd];
            }
            else if (((ckt.Mode & Circuit.Modes.InitJct) != 0) &&
                  ((ckt.Mode & Circuit.Modes.TranOp) != 0) &&
                  ((ckt.Mode & Circuit.Modes.Uic) != 0))
            {
                vds = model.JFETtype * JFETicVDS;
                vgs = model.JFETtype * JFETicVGS;
                vgd = vgs - vds;
            }
            else if (((ckt.Mode & Circuit.Modes.InitJct) != 0) &&
                  (JFEToff == false))
            {
                vgs = -1;
                vgd = -1;
            }
            else if (((ckt.Mode & Circuit.Modes.InitJct) != 0) ||
                  (((ckt.Mode & Circuit.Modes.InitFix) != 0) && (JFEToff)))
            {
                vgs = 0;
                vgd = 0;
            }
            else
            {
                /*
                 *  compute new nonlinear branch voltages 
                 */
                vgs = model.JFETtype *
                    (state.Solution[JFETgateNode] -
                    state.Solution[JFETsourcePrimeNode]);
                vgd = model.JFETtype *
                    (state.Solution[JFETgateNode] -
                    state.Solution[JFETdrainPrimeNode]);

                delvgs = vgs - state.States[0][JFETstate + JFETvgs];
                delvgd = vgd - state.States[0][JFETstate + JFETvgd];
                delvds = delvgs - delvgd;
                cghat = state.States[0][JFETstate + JFETcg] +
                        state.States[0][JFETstate + JFETggd] * delvgd +
                        state.States[0][JFETstate + JFETggs] * delvgs;
                cdhat = state.States[0][JFETstate + JFETcd] +
                        state.States[0][JFETstate + JFETgm] * delvgs +
                        state.States[0][JFETstate + JFETgds] * delvds -
                        state.States[0][JFETstate + JFETggd] * delvgd;
                /*
                 *  limit nonlinear branch voltages 
                 */
                ichk1 = true;
                vgs = pnjlim(vgs, state.States[0][JFETstate + JFETvgs],
                        (JFETtemp * Circuit.CONSTKoverQ), JFETvcrit, ref icheck);
                vgd = pnjlim(vgd, state.States[0][JFETstate + JFETvgd],
                        (JFETtemp * Circuit.CONSTKoverQ), JFETvcrit, ref ichk1);
                if (ichk1)
                {
                    icheck = true;
                }
                vgs = fetlim(vgs, state.States[0][JFETstate + JFETvgs],
                        model.JFETthreshold);
                vgd = fetlim(vgd, state.States[0][JFETstate + JFETvgd],
                        model.JFETthreshold);
            }
            /*
             *   determine dc current and derivatives 
             */
            vds = vgs - vgd;
            if (vgs <= -5 * JFETtemp * Circuit.CONSTKoverQ)
            {
                ggs = -csat / vgs + ckt.Config.Gmin;
                cg = ggs * vgs;
            }
            else
            {
                evgs = Math.Exp(vgs / (JFETtemp * Circuit.CONSTKoverQ));
                ggs = csat * evgs / (JFETtemp * Circuit.CONSTKoverQ) + ckt.Config.Gmin;
                cg = csat * (evgs - 1) + ckt.Config.Gmin * vgs;
            }
            if (vgd <= -5 * (JFETtemp * Circuit.CONSTKoverQ))
            {
                ggd = -csat / vgd + ckt.Config.Gmin;
                cgd = ggd * vgd;
            }
            else
            {
                evgd = Math.Exp(vgd / (JFETtemp * Circuit.CONSTKoverQ));
                ggd = csat * evgd / (JFETtemp * Circuit.CONSTKoverQ) + ckt.Config.Gmin;
                cgd = csat * (evgd - 1) + ckt.Config.Gmin * vgd;
            }
            cg = cg + cgd;

            /* Modification for Sydney University JFET model */
            vto = model.JFETthreshold;
            if (vds >= 0)
            {
                vgst = vgs - vto;
                /*
                 * compute drain current and derivatives for normal mode
                 */
                if (vgst <= 0)
                {
                    /*
                     * normal mode, cutoff region
                     */
                    cdrain = 0;
                    gm = 0;
                    gds = 0;
                }
                else
                {
                    betap = beta * (1 + model.JFETlModulation * vds);
                    Bfac = model.JFETbFac;
                    if (vgst >= vds)
                    {
                        /*
                         * normal mode, linear region
                         */
                        apart = 2 * model.JFETb + 3 * Bfac * (vgst - vds);
                        cpart = vds * (vds * (Bfac * vds - model.JFETb) + vgst * apart);
                        cdrain = betap * cpart;
                        gm = betap * vds * (apart + 3 * Bfac * vgst);
                        gds = betap * (vgst - vds) * apart
                              + beta * model.JFETlModulation * cpart;
                    }
                    else
                    {
                        Bfac = vgst * Bfac;
                        gm = betap * vgst * (2 * model.JFETb + 3 * Bfac);
                        /*
                         * normal mode, saturation region
                         */
                        cpart = vgst * vgst * (model.JFETb + Bfac);
                        cdrain = betap * cpart;
                        gds = model.JFETlModulation * beta * cpart;
                    }
                }
            }
            else
            {
                vgdt = vgd - vto;
                /*
                 * compute drain current and derivatives for inverse mode
                 */
                if (vgdt <= 0)
                {
                    /*
                     * inverse mode, cutoff region
                     */
                    cdrain = 0;
                    gm = 0;
                    gds = 0;
                }
                else
                {
                    betap = beta * (1 - model.JFETlModulation * vds);
                    Bfac = model.JFETbFac;
                    if (vgdt + vds >= 0)
                    {
                        /*
                         * inverse mode, linear region
                         */
                        apart = 2 * model.JFETb + 3 * Bfac * (vgdt + vds);
                        cpart = vds * (-vds * (-Bfac * vds - model.JFETb) + vgdt * apart);
                        cdrain = betap * cpart;
                        gm = betap * vds * (apart + 3 * Bfac * vgdt);
                        gds = betap * (vgdt + vds) * apart
                             - beta * model.JFETlModulation * cpart - gm;
                    }
                    else
                    {
                        Bfac = vgdt * Bfac;
                        gm = -betap * vgdt * (2 * model.JFETb + 3 * Bfac);
                        /*
                         * inverse mode, saturation region
                         */
                        cpart = vgdt * vgdt * (model.JFETb + Bfac);
                        cdrain = -betap * cpart;
                        gds = model.JFETlModulation * beta * cpart - gm;
                    }
                }
            }

            /*
             *   compute equivalent drain current source 
             */
            cd = cdrain - cgd;
            if (((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac | Circuit.Modes.InitSmSig)) != 0) ||
                    (((ckt.Mode & Circuit.Modes.TranOp) != 0) && ((ckt.Mode & Circuit.Modes.Uic) != 0)))
            {
                /* 
                 *    charge storage elements 
                 */
                czgs = JFETtCGS * JFETarea;
                czgd = JFETtCGD * JFETarea;
                twop = JFETtGatePot + JFETtGatePot;
                fcpb2 = JFETcorDepCap * JFETcorDepCap;
                czgsf2 = czgs / model.JFETf2;
                czgdf2 = czgd / model.JFETf2;
                if (vgs < JFETcorDepCap)
                {
                    sarg = Math.Sqrt(1 - vgs / JFETtGatePot);
                    state.States[0][JFETstate + JFETqgs] = twop * czgs * (1 - sarg);
                    capgs = czgs / sarg;
                }
                else
                {
                    state.States[0][JFETstate + JFETqgs] = czgs * JFETf1 +
                            czgsf2 * (model.JFETf3 * (vgs -
                            JFETcorDepCap) + (vgs * vgs - fcpb2) /
                            (twop + twop));
                    capgs = czgsf2 * (model.JFETf3 + vgs / twop);
                }
                if (vgd < JFETcorDepCap)
                {
                    sarg = Math.Sqrt(1 - vgd / JFETtGatePot);
                    state.States[0][JFETstate + JFETqgd] = twop * czgd * (1 - sarg);
                    capgd = czgd / sarg;
                }
                else
                {
                    state.States[0][JFETstate + JFETqgd] = czgd * JFETf1 +
                            czgdf2 * (model.JFETf3 * (vgd -
                            JFETcorDepCap) + (vgd * vgd - fcpb2) /
                            (twop + twop));
                    capgd = czgdf2 * (model.JFETf3 + vgd / twop);
                }
                /*
                 *   store small-signal parameters 
                 */
                if (((ckt.Mode & Circuit.Modes.TranOp) == 0) ||
                        ((ckt.Mode & Circuit.Modes.Uic) == 0))
                {
                    if ((ckt.Mode & Circuit.Modes.InitSmSig) != 0)
                    {
                        state.States[0][JFETstate + JFETqgs] = capgs;
                        state.States[0][JFETstate + JFETqgd] = capgd;
                        return; /*go to 1000*/
                    }
                    /*
                     *   transient analysis 
                     */
                    if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
                    {
                        state.States[1][JFETstate + JFETqgs] =
                                state.States[0][JFETstate + JFETqgs];
                        state.States[1][JFETstate + JFETqgd] =
                                state.States[0][JFETstate + JFETqgd];
                    }
                    var result = ckt.Integration.Integrate(state, JFETstate + JFETqgs, capgs);
                    ggs = ggs + result.Geq;
                    cg = cg + state.States[0][JFETstate + JFETcqgs];
                    result = ckt.Integration.Integrate(state, JFETstate + JFETqgd, capgd);
                    ggd = ggd + result.Geq;
                    cg = cg + state.States[0][JFETstate + JFETcqgd];
                    cd = cd - state.States[0][JFETstate + JFETcqgd];
                    cgd = cgd + state.States[0][JFETstate + JFETcqgd];
                    if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
                    {
                        state.States[1][JFETstate + JFETcqgs] =
                                state.States[0][JFETstate + JFETcqgs];
                        state.States[1][JFETstate + JFETcqgd] =
                                state.States[0][JFETstate + JFETcqgd];
                    }
                }
            }
            /*
             *  check convergence 
             */
            if (((ckt.Mode & Circuit.Modes.InitFix) == 0) || ((ckt.Mode & Circuit.Modes.Uic) == 0))
            {
                if ((icheck)
                        || (Math.Abs(cghat - cg) >= ckt.Config.RelTol *
                            Math.Max(Math.Abs(cghat), Math.Abs(cg)) + ckt.Config.AbsTol) ||
                        (Math.Abs(cdhat - cd) > ckt.Config.RelTol *
                            Math.Max(Math.Abs(cdhat), Math.Abs(cd)) + ckt.Config.AbsTol)
                        )
                {
                    ckt.IsCon = false;
                }
            }
            state.States[0][JFETstate + JFETvgs] = vgs;
            state.States[0][JFETstate + JFETvgd] = vgd;
            state.States[0][JFETstate + JFETcg] = cg;
            state.States[0][JFETstate + JFETcd] = cd;
            state.States[0][JFETstate + JFETcgd] = cgd;
            state.States[0][JFETstate + JFETgm] = gm;
            state.States[0][JFETstate + JFETgds] = gds;
            state.States[0][JFETstate + JFETggs] = ggs;
            state.States[0][JFETstate + JFETggd] = ggd;
            /*
             *    load current vector
             */
            // load:
            ceqgd = model.JFETtype * (cgd - ggd * vgd);
            ceqgs = model.JFETtype * ((cg - cgd) - ggs * vgs);
            cdreq = model.JFETtype * ((cd + cgd) - gds * vds - gm * vgs);
            state.Rhs[JFETgateNode] += (-ceqgs - ceqgd);
            state.Rhs[JFETdrainPrimeNode] +=
                    (-cdreq + ceqgd);
            state.Rhs[JFETsourcePrimeNode] +=
                    (cdreq + ceqgs);
            /*
             *    load y matrix 
             */
            state.Matrix[JFETdrainNode, JFETdrainPrimeNode] += (-gdpr);
            state.Matrix[JFETgateNode, JFETdrainPrimeNode] += (-ggd);
            state.Matrix[JFETgateNode, JFETsourcePrimeNode] += (-ggs);
            state.Matrix[JFETsourceNode, JFETsourcePrimeNode] += (-gspr);
            state.Matrix[JFETdrainPrimeNode, JFETdrainNode] += (-gdpr);
            state.Matrix[JFETdrainPrimeNode, JFETgateNode] += (gm - ggd);
            state.Matrix[JFETdrainPrimeNode, JFETsourcePrimeNode] += (-gds - gm);
            state.Matrix[JFETsourcePrimeNode, JFETgateNode] += (-ggs - gm);
            state.Matrix[JFETsourcePrimeNode, JFETsourceNode] += (-gspr);
            state.Matrix[JFETsourcePrimeNode, JFETdrainPrimeNode] += (-gds);
            state.Matrix[JFETdrainNode, JFETdrainNode] += (gdpr);
            state.Matrix[JFETgateNode, JFETgateNode] += (ggd + ggs);
            state.Matrix[JFETsourceNode, JFETsourceNode] += (gspr);
            state.Matrix[JFETdrainPrimeNode, JFETdrainPrimeNode] += (gdpr + gds + ggd);
            state.Matrix[JFETsourcePrimeNode, JFETsourcePrimeNode] += (gspr + gds + gm + ggs);
        }

        /// <summary>
        /// Load the JFet for AC or PZ analysis
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void AcLoad(Circuit ckt, Model mod, object p = null)
        {
            JfetModel model = mod as JfetModel;

            double gdpr;
            double gspr;
            double gm;
            double gds;
            double ggs;
            double xgs;
            double ggd;
            double xgd;
            CircuitStateComplex state = ckt.StateComplex;
            CircuitState cstate = ckt.State;

            gdpr = model.JFETdrainConduct * JFETarea; // [EDIT] Changed model.JFETdrainResist to model.JFETdrainConduct
            gspr = model.JFETsourceConduct * JFETarea; // [EDIT] Changed model.JFETsourceResist to model.JFETsourceConduct
            gm = cstate.States[0][JFETstate + JFETgm];
            gds = cstate.States[0][JFETstate + JFETgds];
            ggs = cstate.States[0][JFETstate + JFETggs];
            xgs = cstate.States[0][JFETstate + JFETqgs];
            ggd = cstate.States[0][JFETstate + JFETggd];
            xgd = cstate.States[0][JFETstate + JFETqgd];
            state.Matrix[JFETdrainNode, JFETdrainNode] += gdpr;
            state.Matrix[JFETgateNode, JFETgateNode] += ggd + ggs;
            state.Matrix[JFETgateNode, JFETgateNode] += (xgd + xgs) * state.Omega;
            state.Matrix[JFETsourceNode, JFETsourceNode] += gspr;
            state.Matrix[JFETdrainPrimeNode, JFETdrainPrimeNode] += gdpr + gds + ggd;
            state.Matrix[JFETdrainPrimeNode, JFETdrainPrimeNode] += xgd * state.Omega;
            state.Matrix[JFETsourcePrimeNode, JFETsourcePrimeNode] += gspr + gds + gm + ggs;
            state.Matrix[JFETsourcePrimeNode, JFETsourcePrimeNode] += xgs * state.Omega;
            state.Matrix[JFETdrainNode, JFETdrainPrimeNode] -= gdpr;
            state.Matrix[JFETgateNode, JFETdrainPrimeNode] -= ggd;
            state.Matrix[JFETgateNode, JFETdrainPrimeNode] -= xgd * state.Omega;
            state.Matrix[JFETgateNode, JFETsourcePrimeNode] -= ggs;
            state.Matrix[JFETgateNode, JFETsourcePrimeNode] -= xgs * state.Omega;
            state.Matrix[JFETsourceNode, JFETsourcePrimeNode] -= gspr;
            state.Matrix[JFETdrainPrimeNode, JFETdrainNode] -= gdpr;
            state.Matrix[JFETdrainPrimeNode, JFETgateNode] += (-ggd + gm);
            state.Matrix[JFETdrainPrimeNode, JFETgateNode] -= xgd * state.Omega;
            state.Matrix[JFETdrainPrimeNode, JFETsourcePrimeNode] += (-gds - gm);
            state.Matrix[JFETsourcePrimeNode, JFETgateNode] += (-ggs - gm);
            state.Matrix[JFETsourcePrimeNode, JFETgateNode] -= xgs * state.Omega;
            state.Matrix[JFETsourcePrimeNode, JFETsourceNode] -= gspr;
            state.Matrix[JFETsourcePrimeNode, JFETdrainPrimeNode] -= gds;
        }

        /// <summary>
        /// Perform noise analysis
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="mode"></param>
        /// <param name="operation"></param>
        /// <param name="data"></param>
        /// <param name="OnDens"></param>
        /// <param name="p"></param>
        public override void NoiseLoad(Circuit ckt,
            Noise.OutputMode mode,
            Noise.Operation operation,
            NoiseData data,
            ref double OnDens,
            Model mod, object p = null)
        {
            JfetModel model = mod as JfetModel;

            double tmp = 0.0;
            double tempOnoise;
            double tempInoise;
            double[] noizDens = new double[JFETnVar.Length];
            double[] lnNdens = new double[JFETnVar.Length];
            int i;
            CircuitState state = ckt.State;
            NoiseJob curjob = ckt.CurrentJob as NoiseJob;
            string[] JFETnNames = new string[] {       /* Note that we have to keep the order */
	                "_rd",              /* noise due to rd */        /* consistent with the index definitions */
	                "_rs",              /* noise due to rs */        /* in JFETdefs.h */
	                "_id",              /* noise due to id */
	                "_1overf",          /* flicker (1/f) noise */
	                ""                  /* total transistor noise */
                };

            switch (operation)
            {

                case Noise.Operation.N_OPEN:

                    /* see if we have to to produce a summary report */
                    /* if so, name all the noise generators */

                    if (curjob.NStpsSm != 0)
                    {
                        switch (mode)
                        {

                            case Noise.OutputMode.N_DENS:
                                for (i = 0; i < JFETnVar.Length; i++)
                                {
                                    data.namelist.Add(String.Format("onoise_{0}{1}", Name, JFETnNames[i]));
                                }
                                break;

                            case Noise.OutputMode.INT_NOIZ:
                                for (i = 0; i < JFETnVar.Length; i++)
                                {
                                    data.namelist.Add(String.Format("onoise_total_{0}{1}", Name, JFETnNames[i]));
                                    data.namelist.Add(String.Format("inoise_total_{0}{1}", Name, JFETnNames[i]));
                                }
                                break;
                        }
                    }
                    break;

                case Noise.Operation.N_CALC:
                    switch (mode)
                    {

                        case Noise.OutputMode.N_DENS:
                            NevalSrc(ckt, ref noizDens[JFETRDNOIZ], ref lnNdens[JFETRDNOIZ],
                                 Noise.Sources.THERMNOISE, JFETdrainPrimeNode, JFETdrainNode,
                                 model.JFETdrainConduct * JFETarea);

                            NevalSrc(ckt, ref noizDens[JFETRSNOIZ], ref lnNdens[JFETRSNOIZ],
                                 Noise.Sources.THERMNOISE, JFETsourcePrimeNode,
                                 JFETsourceNode, model.JFETsourceConduct * JFETarea);

                            NevalSrc(ckt, ref noizDens[JFETIDNOIZ], ref lnNdens[JFETIDNOIZ],
                                 Noise.Sources.THERMNOISE, JFETdrainPrimeNode,
                                 JFETsourcePrimeNode,
                                 (2.0 / 3.0 * Math.Abs(state.States[0][JFETstate + JFETgm])));

                            NevalSrc(ckt, ref noizDens[JFETFLNOIZ], ref tmp,
                                 Noise.Sources.N_GAIN, JFETdrainPrimeNode,
                                 JFETsourcePrimeNode, (double)0.0);
                            noizDens[JFETFLNOIZ] *= model.JFETfNcoef *
                                 Math.Exp(model.JFETfNexp *
                                 Math.Log(Math.Max(Math.Abs(state.States[0][JFETstate + JFETcd]), Noise.N_MINLOG))) /
                                 data.freq;
                            lnNdens[JFETFLNOIZ] =
                                 Math.Log(Math.Max(noizDens[JFETFLNOIZ], Noise.N_MINLOG));

                            noizDens[JFETTOTNOIZ] = noizDens[JFETRDNOIZ] +
                                             noizDens[JFETRSNOIZ] +
                                             noizDens[JFETIDNOIZ] +
                                             noizDens[JFETFLNOIZ];
                            lnNdens[JFETTOTNOIZ] =
                                 Math.Log(Math.Max(noizDens[JFETTOTNOIZ], Noise.N_MINLOG));

                            OnDens += noizDens[JFETTOTNOIZ];

                            if (data.delFreq == 0.0)
                            {

                                /* if we haven't done any previous integration, we need to */
                                /* initialize our "history" variables                      */

                                for (i = 0; i < JFETnVar.Length; i++)
                                {
                                    JFETnVar[i].LNLSTDENS = lnNdens[i];
                                }

                                /* clear out our integration variables if it's the first pass */

                                if (data.freq == curjob.NstartFreq)
                                {
                                    for (i = 0; i < JFETnVar.Length; i++)
                                    {
                                        JFETnVar[i].OUTNOIZ = 0.0;
                                        JFETnVar[i].INNOIZ = 0.0;
                                    }
                                }
                            }
                            else
                            {   /* data.delFreq != 0.0 (we have to integrate) */
                                for (i = 0; i < JFETnVar.Length; i++)
                                {
                                    if (i != JFETTOTNOIZ)
                                    {
                                        tempOnoise = NoiseIntegration.Integrate(noizDens[i], lnNdens[i],
                                              JFETnVar[i].LNLSTDENS, data);
                                        tempInoise = NoiseIntegration.Integrate(noizDens[i] * data.GainSqInv,
                                              lnNdens[i] + data.lnGainInv,
                                              JFETnVar[i].LNLSTDENS + data.lnGainInv,
                                              data);
                                        JFETnVar[i].LNLSTDENS = lnNdens[i];
                                        data.outNoiz += tempOnoise;
                                        data.inNoise += tempInoise;
                                        if (curjob.NStpsSm != 0)
                                        {
                                            JFETnVar[i].OUTNOIZ += tempOnoise;
                                            JFETnVar[JFETTOTNOIZ].OUTNOIZ += tempOnoise;
                                            JFETnVar[i].INNOIZ += tempInoise;
                                            JFETnVar[JFETTOTNOIZ].INNOIZ += tempInoise;
                                        }
                                    }
                                }
                            }
                            if (data.prtSummary)
                            {
                                for (i = 0; i < JFETnVar.Length; i++)
                                {     /* print a summary report */
                                    data.outpVector[data.outNumber++] = noizDens[i];
                                }
                            }
                            break;

                        case Noise.OutputMode.INT_NOIZ:        /* already calculated, just output */
                            if (curjob.NStpsSm != 0)
                            {
                                for (i = 0; i < JFETnVar.Length; i++)
                                {
                                    data.outpVector[data.outNumber++] = JFETnVar[i].OUTNOIZ;
                                    data.outpVector[data.outNumber++] = JFETnVar[i].INNOIZ;
                                }
                            }    /* if */
                            break;
                    }    /* switch (mode) */
                    break;

                case Noise.Operation.N_CLOSE:
                    return;         /* do nothing, the main calling routine will close */

            }    /* switch (operation) */
        }
    }
}
